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ABSTRACT 


For  purposes  of  studying  the  lateral  heterogeneity  as  well  as  predicting  seismograms  for  this  region,  we  construct  a 
new  3-D  S-velocity  model  by  jointly  inverting  a  variety  of  different  seismic  data.  We  jointly  invert  regional 
waveforms,  surface  wave  group  velocity  measurements,  teleseismic  S  and  SKS  arrival  times,  and  crustal  thickness 
estimates  from  receiver  functions,  refraction  lines,  and  gravity  surveys.  These  data  types  have  complementary 
resolving  power  for  crust  and  mantle  structures,  vertical  and  lateral  variations,  shallow  and  deep  mantle  features, 
local  and  global  structure.  Therefore,  a  joint  inversion  of  these  data  sets  might  help  unravel  the  complexity  of  this 
tectonically  diverse  area.  These  measurements  are  made  from  a  combination  of  MIDSEA,  PASSCAL,  GeoScope, 
Geofon,  GSN,  IDA,  MedNet,  national  networks,  and  local  deployments  throughout  the  study  region  which  extends 
from  the  western  Mediterranean  region  to  the  Hindu  Kush  and  encompasses  northern  Africa,  the  Arabian  peninsula, 
the  Middle  East,  and  part  of  the  Atlantic  Ocean  for  reference.  The  Moho  depth  result  is  broadly  consistent  with 
CRUST2.0,  except  in  mid-northern  Africa,  where  the  crust  from  our  joint  inversion  is  about  5  km  thinner.  Fast 
velocity  anomalies  are  found  beneath  the  West  African  Craton,  the  Hellenic  trench,  the  Apennines,  the  East 
European  Platform,  and  the  Arabian  Platform  at  a  75-150  km  depth,  whereas  low- velocity  anomalies  are  located 
along  the  plate  boundaries  such  as  the  mid- Atlantic  ridge,  Afar,  the  Anatolian  Plateau,  Iran,  Afghanistan,  western 
Mediterranean  Sea,  and  the  Red  Sea.  Based  on  the  assumption  that  P-wave  velocity  anomalies  in  m/s  scale  are  very 
close  to  S-wave  velocity  anomalies  in  m/s  scale  in  our  study  area,  we  convert  our  S-velocity  model  to  a  P- velocity 
model  with  use  of  teleseismic  P  arrival  times.  Finally,  our  model  is  validated  by  performing  travel-time  predictions 
with  a  dataset  of  ground  truth  events.  Our  model  generally  produces  better  travel-time  predictions  than  the  iasp91 
model.  We  expect  that  better  travel-time  predictions  may  be  achieved  by  replacing  the  poorly  resolved  crust  in  our 
model  with  CRUST2.0  and  refining  the  P- velocity  model. 
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OBJECTIVES 

Our  primary  objective  is  developing  a  new  3-D  ^-velocity  model  for  the  Middle  East  and  Mediterranean  region, 
including  North  Africa,  southern  Europe,  and  Arabia  that 

1)  is  resolved  in  aseismic  regions; 

2)  is  resolved  throughout  the  upper  mantle  (to  660  km); 

3)  resolves  laterally  varying  crustal  thickness; 

4)  contains  laterally  varying  vertical  velocity  gradients; 

5)  is  simultaneously  compatible  with  multiple  data  sets; 

6)  utilizes  several  recent,  unique  waveform  data  sets;  and 

7)  includes  uncertainties  of  the  model  parameters. 

These  features  would  increase  the  model’s  ability  to  predict  and  calibrate  regional  travel  times  and  waveforms, 
thereby  providing  improved  event  locations,  focal  mechanisms,  and  other  event  discriminants. 

Secondly,  we  aim  to  convert  the  3-D  S-velocity  model  to  a  3-D  P-velocity  model,  using  teleseismic  P-arrival  times. 

Our  third  objective  is  to  test  both  the  P- and  S-wave  models’  ability  to  predict  regional  P  and  S  travel  times,  deflect 
wave  paths  and  deform  waveforms,  and  assess  their  effects  first  on  the  studied  seismograms  (travel  times  and 
amplitudes)  and  subsequently  on  the  3-D  models  derived  from  these  data. 

This  report  covers  all  the  above-mentioned  objectives. 


Figure  1.  Topographic  map  of  the  study  region.  The  pink  line  is  the  NUVEL1-A  (DeMets  et  al.,  1990) 
representation  of  the  Eurasia-Africa-Arabia  plate  boundary. 

The  study  region  is  centered  around  the  Africa- Arabia-Eurasia  triple  junction  (Figure  1),  and  extends  west  to  the 
Africa-Eurasia-North  America  junction  at  the  Azores,  just  off  the  map,  and  east  to  the  Arabia-Eurasia-Indian  Plate 
junction.  The  NUVEL1-A  (DeMets  et  al.,  1990)  representation  of  these  plate  boundaries  is  shown  by  the  pink  line  in 
Figure  1.  The  interaction  of  these  five  major  tectonic  plates  with  each  other  and  with  several  microplates  within  an 
area  of  one  quarter  of  the  Earth’s  circumference  yields  this  region  rich  with  tectonic  complexity.  We  plan  to  capture 
various  renditions  of  this  structurally  complicated  part  of  the  world  in  one  S-velocity  model  through  the  joint 
inversion  of  several  different  types  of  seismic  data  simultaneously;  the  new  model  will  refine  our  understanding  of 
the  structure  and  tectonics  in  this  region  of  the  Earth.  The  data  types  we  combine  are  constraints  from  independent 
studies  on  the  depth  to  the  Moho,  fundamental-mode  Rayleigh-wave  group  velocity  measurements,  waveform  fits  of 
regional  S  and  Rayleigh  waves,  and  arrival  times  of  teleseismic  S  and  SKS  waves.  We  convert  the  resulting  S- 
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velocity  model  to  a  ^-velocity  model  with  over  2.9  million  teleseismic  P  arrival  times.  Then,  we  validate  our  model 
by  performing  travel-time  prediction  with  a  dataset  of  ground  truth  event. 

RESEARCH  ACCOMPLISHED 


(c)  Group  velocities  (d)  Moho  depth  constraints 


Figure  2.  Datasets  used  for  the  joint  inversion,  (a)  Ray  path  coverage  for  regional  waveform  fits.  Stations  are 
illustrated  as  red  triangles,  and  events  as  yellow  circles,  (b)  Events  and  stations  used  for  teleseismic 
S  and  SKS  arrival  time  estimation.  Cyan  circles  and  blue  triangles  represent  events  and  stations 
used  for  relative  delay  time  estimation  with  MCCC,  respectively.  Yellow  circles  and  red  triangles 
are  events  and  stations  from  the  reprocessed  ISC  catalogue,  (c)  Great-circle  wave  paths  for  45  s 
period  Rayleigh  waves.  Stations  are  illustrated  as  red  triangles,  and  events  as  yellow  circles,  (d)  Map 
of  the  Moho  depth  distribution  acquired  from  literatures.  Artificial  point  constraints  of  10  km  depth 
are  placed  to  the  Atlantic  and  Indian  Oceans  where  measurements  are  absent. 


Regional  Waveform  Fits 

We  fit  over  5,600  available  waveforms  from  the  Lawrence  Livermore  National  Laboratory  (LLNL)  database  and 
MIDSEA  dataset  (Marone  et  al.,  2004)  which  sample  the  Mediterranean  region,  North  Africa,  the  Middle  East, 
Pakistan,  and  Afghanistan  using  the  non-linear  inversion  procedure  employed  by  partitioned  waveform  inversion 
studies  (Nolet,  1990;  Van  der  Lee  and  Nolet,  1997).  We  utilize  events  with  magnitude  larger  than  4.0  and 
seismograms  with  epicentral  distance  from  5°  to  50°.  The  great-circle  wave  paths  for  these  seismograms  are  shown 
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in  Figure  2(a).  We  have  estimated  path-averaged  S- velocity  structures  for  these  paths  using  the  same  starting  S- 
velocity  model  but  a  different  crustal  thickness  (in  5  km  increments)  for  each  path,  based  on  a  priori  reported 
estimates  (Marone  et  ah,  2004).  The  starting  model  often  predicts  significant  phase  differences  relative  to  the  data 
for  both  the  S-  and  Rayleigh  waves.  The  inversion  procedure  estimates  the  perturbations  to  the  starting  model  by 
non-linear  optimization  (Nolet,  1990;  van  der  Lee  and  Nolet,  1997).  For  the  computation  of  synthetic  waveforms  we 
also  need  to  know  the  source  mechanisms  of  the  events.  We  typically  use  the  centroid  moment  tensors  (CMT)  from 
the  Global  CMT  project,  and  hypocenters  from  Engdahl’s  reprocessed  International  Seismological  Centre  (ISC) 
database  (Engdahl  et  al.,  1998).  The  frequency  content  of  the  data  and  synthetic  is  different  for  each  fit  and 
generally  falls  within  the  band  0.006  to  0.1  Hz. 

Teleseismic  S  and  SKS  Arrival  Times 

We  obtained  S  and  SKS  phase  arrival  time  data  from  two  different  sources.  Both  arrival  times  are  adjusted  for 
topography  and  Earth’s  ellipticity  before  inversion.  First,  we  used  high-quality  relative  arrival  times  of  teleseismic  S 
and  SKS  waves  (Benoit  et  al.,  2006;  Park  et  al.,  2007;  Schmid  et  al.,  2004)  which  are  measured  on  broadband 
seismograms  with  use  of  multi-channel  cross-correlation  (MCCC;  VanDecar  and  Crosson,  1990).  Benoit  et  al. 
(2006),  Park  et  al.  (2007),  and  Schmid  et  al.  (2004)  used  seismograms  recorded  in  the  Mediterranean,  Ethiopia,  amd 
Saudi  Arabia,  respectively.  Moreover,  we  measured  additional  relative  delay  times  at  Turkey  and  central  Asia.  The 
number  of  S  phase  relative  arrival  times  is  over  5900  with  epicentral  distances  of  30°-90°  and  the  number  of  SKS 
phases  is  over  1,400  with  distances  of  87°-140°;  over  73,00  in  total.  We  set  0.5  s  as  uncertainty  for  weights  during 
the  inversion.  Second,  we  obtained  over  223,000  S  phase  arrival  time  data  from  the  reprocessed  ISC  database 
(Engdahl  et  al.,  1998)  from  1964  to  2007.  The  epicentral  distance  range  is  20°-80°.  We  set  lower  uncertainty  of  2.5  s 
for  the  data  because  lots  of  outliers  and  systematic  errors  are  found  in  this  data  set  by  Rohm  et  al.  (2000).  We 
illustrate  location  of  stations  and  events  for  each  of  the  two  types  of  arrival  time  data  in  Figure  2(b). 

Surface  Wave  Group  Velocities 

We  have  measured  over  105,000  group  velocities  of  fundamental-mode  Rayleigh  waves  recorded  at  MIDSEA  and 
other  stations  in  the  region  and  used  them  to  update  previous  group  velocity  maps  (e.g.,  Pasyanos,  2005).  Figure  2(c) 
shows  the  great-circle  wave  paths  for  which  we  included  fundamental-mode  Rayleigh  wave  group  velocities  in  our 
joint  inversion.  Compared  to  other  data  types  we  include,  lateral  coverage  is  best  for  this  group  velocity  data  set, 
though  vertical  coverage  is  provided  mainly  by  the  teleseismic  arrival  times  and  regional  S  and  Rayleigh  waveforms. 
The  period  for  group  velocities  ranges  from  7  to  100  sec. 

Crustal  Thickness  Constraints 

While  the  fundamental-mode  Rayleigh- wave  group  velocities  and  regional  S  and  Rayleigh  waveforms  have 
significant  sensitivity  to  Moho  depth,  they  cannot  uniquely  resolve  Moho  depth.  We  therefore  include  independent 
estimates  of  crustal  thickness  as  point  constraints  in  the  joint  inversion  and  thus  compile  such  measurements  from  a 
large  number  of  published  studies.  A  partial  list  of  these  studies  is  provided  in  Marone  et  al.  (2003)  and  part  of  these 
points  stem  from  the  database  used  for  CRUST5.1  (Mooney  et  al.,  1998).  We  have  added  new  estimates  of  Moho 
depth  from  more  recent  studies  and  interactively  resolved  or  removed  conflicting  data  and  outliers  from  the  data  set. 
Over  4,700  Moho  depth  constraints  are  mapped  in  Figure  2(d).  The  majority  of  this  data  set  is  from  receiver  function 
studies.  Remaining  points  are  from  active-source  studies  and  some  are  from  gravity  surveys.  The  active-source 
constraints  were  assigned  relatively  low  errors  and  the  Moho-depth  estimates  from  gravity  were  assigned  the  largest 
errors.  For  the  oceans  we  use  a  constraint  of  10  km  for  Moho  depth,  but  only  for  points  also  covered  by  data  from 
our  other  data  sets. 

Inversion  Results 

We  have  developed  software  that  handles  the  joint  inversion  of  constraints  from  regional  waveform  fits,  crustal 
estimates,  group  velocities,  and  teleseismic  arrival  times.  The  resolving  power  of  the  combined  data  is  superior  to 
that  of  each  of  the  data  sets  alone.  The  joint  inversion  reduces  the  variance  in  the  data  sets  by  38%  for  the 
teleseismic  delay  times,  55%  for  the  Rayleigh- wave  group  velocities,  90%  for  the  Moho  point  constraints,  and  88% 
for  the  regional  waveform  fits. 
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Figure  3  shows  the  joint  inversion  results  for  Moho  depth  and  velocity  perturbations.  The  Moho  map  shows  a  good 
resemblance  with  the  point  constraints  in  Figure  2(d),  but  is  much  smoother  because  the  regional  waveforms  and 
group  velocity  data,  as  well  as  regularization  constraints  in  the  joint  inversion  provide  a  smooth  interpolation 
between  the  point  constraints.  The  map  is  also  broadly  consistent  with  CRUST2.0,  except  in  northern  Africa,  where 
the  crust  from  our  joint  inversion  is  about  5  km  thinner.  Resolution  tests  indicate  that  our  combined  data  sets  can 
resolve  the  crustal  thicknesses  of  CRUST2.0. 

Fast-velocity  anomalies  are  found  beneath  the  West  African  Craton,  the  Hellenic  trench,  the  Apennines,  the  East 
European  Platform,  and  the  Arabian  Platform  at  75-150  km  depth,  whereas  low-velocity  anomalies  are  located 
along  the  plate  boundaries  such  as  the  mid- Atlantic  ridge,  Afar,  the  Anatolian  Plateau,  Iran,  Afghanistan,  western 
Mediterranean  Sea,  and  the  Red  Sea.  Fast-velocity  anomalies  are  spread  beneath  the  western  Mediterranean  Sea, 
Dinarides,  the  Aegean  Sea,  Turkey,  and  the  Black  Sea  at  a  500-km  depth,  which  may  indicate  subducting  slabs  from 
the  Apennine,  the  Calarian  arc,  and  the  Hellenic  arc. 


(c)  300  km  depth  (d)  500  km  depth 

Figure  3.  Moho  and  velocity  perturbations  resulting  from  our  joint  inversion.  Moho  map  (a)  and  velocity 
perturbations  at  100  km  (b),  300  km  (c),  and  500  km  depth  (d)  are  illustrated. 
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Conversion  of  N-Velocity  Model  to  P-Velocity  Model 

One  way  to  obtain  a  ^-velocity  model  is  to  perform  P-wave  arrival  time  inversion  using  the  equation, 

Gama=da,  (1) 

where  Ga  is  the  sensitivity  kernel  matrix  for  P-wave  arrival  times,  ma  is  the  3-D  ^-velocity  model  vector,  and 

dais  the  teleseismic  P-wave  delay  vector  which  is  obtained  by  subtracting  predicted  arrival  times  through  a 

reference  model  from  observed  arrival  times.  However,  it  is  difficult  to  resolve  shallow  structure  with  teleseismic  P- 
wave  arrival  time  data  because  of  poor  resolution  for  shallow  depth.  Therefore,  we  have  to  find  another  way  to  get 
P- velocity  model  with  P- wave  arrival  time  data.  It  would  be  very  promising  in  the  inversion  for  P- velocity  model  if 
we  can  adopt  the  S-velocity  model  as  a  preliminary  model,  because  the  S-velocity  model  has  good  resolution  for 
both  shallow  and  deep  structure  due  to  our  use  of  regional  waveform  fits,  teleseismic  arrival  times,  Rayleigh-wave 
group  velocities,  and  Moho  depth  constraints. 


We  assume  that  P-wave  velocity  anomalies  in  m/s  scale  are  very  close  to  the  S-wave  velocity  anomalies  in  m/s  scale 
in  our  study  area  if  temperature  dominates  the  for  cause  of  the  velocity  perturbations  based  on  results  in  et  al.  (2004) 
in  which  similarity  between  P-  and  S-  velocity  anomalies  is  drawn  with  use  of  delay  time  ratio  (see  Figure  4a). 
Therefore,  it  is  reasonable  to  adopt  our  S- velocity  model  as  a  starting  model  in  performing  a  P- wave  arrival  time 
inversion  (Schmid  et  al.,  2008).  By  multiplication  of  the  S-velocity  model  and  EAvave  sensitivity  kernels,  we  can  get 

predicted  P- wave  delays  d 6^a  by  as  in 

G„nv  =  dfi_a ,  (2) 

where  is  the  3-D  ^-velocity  model  in  m/s  scale,  and  d^a  is  predicted  P  delays  by  .  Then,  teleseismic  P- 
wave  delays  d  a  can  be  divided  as 

d«  =  +  dre, ,  (3) 

where  dres  is  the  residual  P  delay  vector  which  means  the  difference  between  da  and  d^a  .  Therefore,  Eq.  (1) 
can  be  rewritten  as 

Ga(mfi+mnew)  =  d^a+dres-  (4) 

Eq.  (4)  is  a  linear  equation,  so  we  extract  Eq.  (2)  from  Eq.  (4)  to  perform  P- wave  arrival  time  inversion  for  mnew 
using  following  equation, 

=  dres .  (5) 

A  P-  velocity  model  is  finally  obtained  by  adding  mwew  to  the  S-velocity  model, 


m«  =  +mnew.  (6) 

The  resulting  P- velocity  model  with  teleseismic  P  arrival  time  data  in  Figure  4(b)  is  shown  at  various  depths  in 
Figure  5.  The  results  are  similar  to  S-velocity  perturbation  results  in  Figure  3,  which  means  P-velocity  variations 
from  the  S-velocity  model  are  small.  But  higher-velocity  anomalies  are  found  in  the  Alps,  Hellenides,  Turkey,  and 
the  Hindu  Kush  than  in  the  S-velocity  model. 


Model  Evaluation 


Distributions  of  the  travel-time  residuals  with  respect  to  iasp91  and  our  joint  inversion  model  (Joint  Model)  are 
shown  in  histogram  form  in  Figure  6  for  stations  AJM,  ELL,  FRU,  MAIO,  NIL,  PTO,  QUE,  and  SVE.  We  report 
both  the  L2  norm  (mean,  standard  deviation,  STD  and  variance  reduction,  VR)  and  LI  norm  statistics  (median, 
scaled  median  absolute  deviation,  SMAD  and  SMAD  reduction,  SMADR).  The  SMAD  provides  an  estimate  of  the 
spread  of  values  which  is  less  sensitive  to  outliers  and  so  is  often  more  appropriate  in  the  presence  of  non-Gaussian 
errors.  At  most  stations,  SMADR  is  positive  except  for  ELL,  MAIO  and  QUE. 
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(a)  (b) 

Figure  4.  (a)  S  delays  versus  P  delays  (Schmid  et  al.,  2004).  (b)  Events  and  stations  used  for  teleseismic  P 
arrival  time  estimation.  Cyan  circles  and  blue  triangles  represent  events  and  stations  used  for 
relative  delay  time  estimation  with  MCCC,  respectively.  Yellow  circles  and  red  triangles  mean 
events  and  stations  from  the  reprocessed  ISC  catalogue. 


(a)  100  km  depth  (b)  200  km  depth 

Figure  5.  P-velocity  perturbations  at  (a)  100  and  (b)  200  km  depth. 

Visually  we  can  evaluate  how  well  the  joint  inversion  model  (Joint  Model)  predicts  the  data  geographically  by 
computing  3D  travel-time  correction  surfaces  for  some  stations.  To  compute  such  model-based  correction  surfaces 
we  subtract  the  iasp91 -predicted  time  from  the  Joint  Model-predicted  time  using  3D  finite-difference  algorithm. 
Example  surfaces  are  shown  in  Figure  6  for  the  stations  for  a  source  depth  of  10  km;  we  choose  this  depth  as  it 
represents  the  average  focal  depth  of  the  GT25  events  (epicenter  mislocation  of  25  km  or  less)  selected.  Color-coded 
residuals  (data-iasp91)  from  events  of  focal  depth  0  to  20  km  are  also  plotted  on  top  of  the  surfaces  in  Figure  6  to 
provide  a  visual  assessment  of  the  geographic  agreement  between  Joint  Model  prediction  and  the  data. 
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Figure  6.  Histograms  of  travel-time  residuals  for  both  models  at  stations  AJM,  ELL,  FRU,  MAIO,  NIL,  PTO, 
QUE,  and  SVE  along  with  the  travel-time  residual  surfaces  (Joint  Model-iasp91). 
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CONCLUSIONS  AND  RECOMMENDATIONS 


Through  a  joint  inversion  of  teleseismic  S-wave  arrival  time  delays,  waveform  fits  of  regional  S-  and  Rayleigh 
waves,  group  velocity  measurements  of  fundamental-mode  Rayleigh  waves,  and  independent  constraints  on  Moho 
depth,  we  have  achieved  considerable  variance  reduction  in  each  data  set  simultaneously.  The  new  3-D  S-velocity 
model  is  converted  to  a  P- velocity  model  with  teleseismic  P  arrival  times  and  the  ^-velocity  model  is  used  for 
travel-time  prediction.  Our  model  generally  produces  better  prediction  than  the  iasp91  ID  model.  SMADR  is  not 
reduced  exceptionally  at  stations  ELL,  MAIO  and  QUE,  but  mean  and  median  are  closer  to  zero  at  station  ELL  and 
QUE  in  Joint  Model  residuals,  which  means  our  model  succeeded  in  finding  average  velocity  perturbations.  Better 
refined  results  could  be  achieved  by  substituting  our  poorly  resolved  crust  with  CRUST2.0  and  refining  the  P- 
velocity  model. 
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